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it can be used in conjunction with norm-conserving, non-local pseudopotentials. This extension 
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on several finite and periodic systems, finding very good agreement with established methods and 
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I. INTRODUCTION 

The experimental technique of nuclear magnetic reso- 
nance (NMR) is a powerful tool to determine the struc- 
ture of molecules, liquids, and periodic systems. It is thus 
not surprising that, since its discovery in 1938, NMR 
has evolved into one of the most widely used methods 
in structural chemistry. ^'^ Unfortunately, one caveat of 
this successful method is that there is no basic, generally 
valid "recipe" that allows a unique determination of the 
structure given a measured spectrum. As a result, for 
more complex systems the mapping between structure 
and measured spectrum can be ambiguous. 

It had been realized early on that ab-initio calcula- 
tions could resolve some of these ambiguities and thus 
greatly aid in determining structures from experimen- 
tal NMR spectra. For finite systems such as simple 
molecules, appropriate methods were first developed in 
the quantum-chemistry community.^ While highly accu- 
rate, these methods by construction were unable to cal- 
culate NMR shifts of periodic systems, which is impor- 
tant for the increasingly popular solid-state NMR spec- 
troscopy. The underlying physical limitation is due to 
the fact that the description of any constant external 
magnetic field requires a non-periodic vector potential. 
Possible approaches to combine such a non-periodic vec- 
tor potential with the periodic potential of crystals were 
found only within the last decade.'*"* All of these ap- 
proaches have in common that they treat the external 
magnetic field in terms of the linear-response it causes 
to the system under consideration. Although these ap- 
proaches are accurate and successful, the required linear- 
response framework makes them fairly complex and dif- 
ficult to implement. 

Recently, a fundamentally different approach for the 
calculation of ab-initio NMR shifts has been developed 
by some of us.^ In our converse approach we circumvent 
the need for a linear-response framework in that we re- 
late the shifts to the macroscopic magnetization induced 
by magnetic point dipoles placed at the nuclear sites of 
interest. The converse approach has the advantage of be- 



ing conceptually much simpler than other standard ap- 
proaches and it also allows us to calculate the NMR shifts 
of systems with several hundred atoms. Our converse 
approach has already successfully been applied to sim- 
ple molecules, crystals, liquids, and extended poly cyclic 
aromatic hydrocarbons. ^'*° 

While the converse method can be directly imple- 
mented in an all-electron first-principles computer code, 
an implementation into a pseudopotential code is more 
complicated due to nonlocal projectors usually used in 
the Kleinman-Bylander separable form.*^ In this pa- 
per we present a mathematical extension to the con- 
verse formalism such that it can be used in conjunction 
with norm-conserving, non-local pseudopotential. This 
extension permits the efficient ab-initio calculation of 
NMR chemical shifts for elements other than hydrogen 
within the convenience of a plane-wave pseudopotential 
approach. 

The paper is organized in the following way: In Sec. II 
we first review the converse-NMR method and its rela- 
tion to the orbital magnetization at an all-electron level. 
Then, we apply the gauge including projector augmented 
wave (GIPAW) transformation to derive an expression 
for the orbital magnetization in the context of norm- 
conserving pseudopotentials. We discuss aspects of the 
implementation of the converse method in Sec. III. In 
this section we also show results of several convergence 
tests that we have performed. To validate our approach, 
we apply our converse approach to molecules and solids 
and the results are collected in Sec. IV. In Sec. V we 
discuss the main advantages of the converse method. Fi- 
nally, we summarize and conclude in Sec. VI. The GI- 
PAW transformation and several details of the mathe- 
matical formalism — which would be distracting in the 
main text — are presented in Appendices A and B. 

II. THEORY 

The converse method for calculating the NMR chemi- 
cal shielding has been introduced in Ref. [9] and can be 
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summarized as follows: 



dm,., 



(1) 



Thus, in the converse method the chemical shielding ten- 
sor CTs, Q/3 is obtained from the derivative of the orbital 
magnetization M with respect to a magnetic point dipole 
rris, placed at the site of atom s. Sap is the Kronecker 
delta and D, is the volume of the simulation cell. In other 
words, instead of applying a constant magnetic field to 
an infinite periodic system and calculating the induced 
field at all equivalent s nuclei, wc apply an infinite array 
of magnetic dipolcs to all equivalent sites s, and calcu- 
late the change in magnetization. Since the perturbation 
is now periodic, the original problem of the non-periodic 
vector potential has been circumvented. 

In practice, the derivative in Eq. (I) is calculated as a 
finite difference of the orbital magnetization in presence 
of a small magnetic point dipole mg. Since M vanishes 
for rris = and is an odd function of nis because of time- 
reversal symmetry (for a non-magnetic system in absence 
of spin-orbit interaction), it is sufficient to perform three 
calculations of M(mse), where e are cartesian unit vec- 
tors. 

Within density-functional theory (DFT) the all- 
electron Hamiltonian is, in atomic units: 



where w„k are the Bloch wavefunctions, 
— e^'^^'' H e*'"', e„k are its eigenvalues and ep 
is the Fermi level. The k-derivative of the Bloch wave- 
functions can be evaluated as a covariant derivative, 
or by k • p perturbation theory. 

Equations (2) and (5) arc adequate to evaluate the 
NMR shielding tensor Eq. (1) in the context of an all- 
electron method (such as FLAPW, or local-basis meth- 
ods), in which the interaction between core and valence 
electrons is treated explicitly. However, in a pseudopo- 
tential framework, where the effect of the core elec- 
trons has been replaced by a smooth effective potential, 
Eqs. (2) and (5) arc not sufficient to evaluate the NMR 
shielding tensor. One reason for this is that the valence 
wave functions have been replaced by smoother pseudo 
wave functions which deviate significantly from the all- 
electron ones in the core region. 

In the following sections, we derive the formulas needed 
to calculate the converse NMR shielding tensor in Eq. (1), 
in the context of the pscudopotcntial method. Our 
derivation is based on the GIPAW transformation^ (see 
also Appendix A), that allows one to reconstruct all- 
electron wave functions from smooth pseudopotential 
wave functions. For the sake of simplicity, we assume 
all GIPAW projectors to be norm-conserving. 



+ Vks(r) 



where 



A«(r) = 



ms X (r - r^) 



'-rJ3 



(2) 



(3) 



is the vector potential corresponding to a magnetic dipole 
nis centered at the atom s coordinate r^.^^ We neglect 
any explicit dependence of the exchange-correlation func- 
tional on the current density. In practice, spin-current 
density- functional theory calculations have shown to pro- 
duce negligible corrections to the orbital magnetiza- 



tion 
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In finite systems (i.e. molecules), the orbital magneti- 
zation can be easily evaluated via the velocity operator 
v = -i[r,H]: 



M = ^^(V'„|rxv|^„ 



(4) 



where l^^^) are molecular orbitals, spanning the occupied 
manifold. For periodic systems, the situation is more 
complicated due to itinerant surface currents and to the 
incompatibility of the position operator r with periodic 
boundary conditions. It has been recently shown^^"^^ 
that the orbital magnetization in a periodic system is 
given by: 



M = -^Im^/(e„k) {dkUnk\ 

nk 



X {Hu + e„k - 2eF) \dkUnu) (5) 



A. The converse-NMR GIPAW hamiltonian 

In this section we derive the pseudopotential GI- 
PAW hamiltonian corresponding to the all-electron (AE) 
hamiltonian in Eq. (2). For reasons that will be clear in 
the next section, we include an external uniform mag- 
netic field B in addition to the magnetic field generated 
by the point dipole uig. For the sake of simplicity, we 
carry out the derivation for an isolated system (i.e. a 
molecule) in the symmetric gauge A(r) = (f/2)B x r. 
The generalization to periodic systems is then performed 
at the end. 

We start with the all-electron hamiltonian: 



p+i(A(r)+A,(r)) 



+ V{r) (6) 



We now decompose Eq. (6) in powers of A^ as T^ae = 



n 



AE 



n 



(si) 
AE 



^AE^ , where 



« = ^(p+^A(r))%nr) 



(7) 



n 



{si) 
AE 

,(s2) 



- (p • A,(r) + A,(r) . p + -A(r) • A,(r))(8) 



We can neglect 1-1^^^ in all calculations, since nis is a 
small perturbation to the electronic structure. 

We then apply the GIPAW transformation Eq. (A4) 
to the two remaining terms, Eqs. (7) and (8), and we ex- 
pand the results up to first order in the magnetic field. ""^ 
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At zeroth order in the external magnetic field B, the 

GIPAW transformation of "Hae' ^^'^ '^ae yields the GI- 
PAW hamiltonian: 



IPAW 



(sO,0) 
GIPAW 



"•GIPAW "-GIPAW 



(9) 
(10) 



^(sl,0) 
"•GIPAW 



R 

i[p.A,(r) + A,(r).p]+5^ifr(ll) 



R 



where 14oc(r) is the local Kohn-Sham potential and 
is the non-local pseudopotential in the separarable 
Kleinmann-Bylander (KB) form: 



(12) 



Similar to V^^, the term has the form of a non-local 
operator 

= (13) 
nm 

kR,nm = (0R,„|p • As(r) + As(r) • p|(/)R,„) - 

(^R,„|p- A,(r)-HAs(r) -pI^r,™) (14) 

The index R runs over all atoms in the system, and the 
indexes n and m, individually run over all projectors as- 
sociated with atom R. For a definition of |<^R,n) and 
|</'R,n) see Appendix A. Note that the set of GIPAW pro- 
jectors Ipr) need not be the same as the KB projectors 
|/3r). For instance, in the case of norm-conserving pseu- 
dopotentials, one KB projector per non-local channel is 
usually constructed. Conversely, two GIPAW projectors 
for each angular momentum channel are usually needed. 

Equation (9) is the Hamiltonian to be implemented in 
order to apply a point magnetic dipole to the system. 
The first term of ^qjpaw ^® applied to a wave func- 
tion in real space or in reciprocal space. The second term 
acts on the wave functions like an extra non-local term 
and requires very little change to the existing framework 
that applies the non-local potential. 

At the first order in the magnetic field, the GIPAW 
transformation yields two terms: 



"GIPAW 



"GIPAW 



lB.(L + $:Rxi[r,yNL]) (15) 

R 

^B.(rxA,(r) + 5^E^i^ + 

R, 

+ ^Rxi[r,i^NL] 



R 



(16) 



where L = r x p and E^'" is the non-local operator 

Er'^ = ^\Pn,n)en,nm{Pn,m\ (17) 

nm 

eR,nm = (0R,n|(r - R) X As{r)\(f)R,m) - 

(0R,„|(r-R) X A,(r)|^R,„) (18) 



The two equations above will be used in the next section, 
in conjunction with the Hellmann-Feynman theorem, to 
derive the GIPAW form of the orbital magnetization. 



The orbital magnetization in the GIPAW 
formalism 



The orbital magnetization for a non spin-polarized sys- 
tem is formally given by the Hellmann-Feynman theorem 
as 



(19) 



In the GIPAW formalism this expectation value can 
be expressed in terms of the GIPAW Hamiltonian and 
pseudo wave functions 

9'Hgipaw 



M = - 



dB 



(20) 



B=0 



By using the results of the previous section we find: 

occ 

M = -2^ E + ^ >^ w + E^R + 

n R 

+ z^Rx [l/^L^ii:|[L,r] U„) (21) 



R 

Note that in the expression above, ipn arc the eigcnstates 
of the GIPAW Hamiltonian in absence of any external 
magnetic fields. 

While the formula (21) for M can directly be applied 
to atoms and molecules, it is ill-defined in the context of 
periodic systems, owing to the presence of the position 
operator — explicitly as in r, but also implicitly as in L. 
This problem can be remedied by applying the modern 
theory of orbital magnetization. The goal is thus to 
reformulate Eq. (21) in terms of (r x vqipaw)- We can 
calculate this operator as 

r X VGiPAw = r X T[r, ■Hgipaw]b=o = 
^ 

= L + rxA,{r) + iJ2rx[Vi^ + K^'',r] (22) 

R 

Replacing L in Eq. (21) by the corresponding expression 
calculated from Eq. (22), and regrouping the terms, we 
obtain the central result of this paper: 



M 



M 



bare 



■ Mnl + M 



para 



■Mdi, 



Mba 



(r X VGIPAW ) 



M 



NL 



Mdi, 




(R-r)x [(R-r),y^ 



NLl 



NLl 
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TABLE I: Chemical shielding a in ppm of a hydrogen atom 
in a water molecule as a function of the kinetic-energy cutoiT 
£kin (in units of Rydberg). 



-Ekin 


a 


i'kin 


a 


(Ry) 


(ppm) 


(Ry) 


(ppm) 


30 


31.0009 


70 


31.1177 


40 


31.0595 


80 


31.1301 


50 


31.0637 


90 


31.1228 


60 


31.0832 


100 


31.1119 



where (. . .) stands for X^n™ (V'nl • • • IVAi)- The naming of 
the various terms are in analogy to Ref. [6]. The set of 
equations (23)-(26) are now valid both in isolated and 
periodic systems, as shown in detail in Appendix B. 



III. IMPLEMENTATION AND 
COMPUTATIONAL DETAILS 

We have implemented the converse NMR method and 
its GIPAW transformation into PWSCF, which is part of 
the Quantum-Espresso package.^° 

In principle, the calculation of the NMR shielding is 
performed the following way: The vector potential cor- 
responding to the microscopic dipole is included in the 
Hamiltonian and the Kohn-Sham equation is solved self- 
consistently under that Hamiltonian. In practice, how- 
ever, in a first step one can equally as well find the ground 
state of the unperturbed system. Based on this ground 
state, in the second step one can then introduce the 
dipole perturbation and reconverge to the new ground 
state. Note that the reconvergence of the small dipole 
perturbation is usually very fast and only a small number 
of SCF steps is necessary in addition to the ground state 
calculation. In fact, tests have shown that a reconver- 
gence is not even necessary — diagonalizing the perturbed 
Hamiltonian only once with the unperturbed wave func- 
tions gives results for the NMR shielding within 0.01 ppm 
of the fully converged solution. This yields a huge cal- 
culational benefit for large systems, where wc calculate 
the unperturbed ground state once and then, based on 
the converged ground-state wave functions, calculate all 
shieldings of interest by non-SCF calculations. 

In order to study the convergence of the NMR chemical 
shielding with several parameters, we performed simple 
tests on a water molecule in the gas phase. The molecule 
was relaxed in a box of 30 Bohr; for all our calculations 
wc^ uscxl Troullicr-Martin norm-conserving psciidopoten- 
tials^^ and a PBE exchange-correlation functional. 

First, we tested the convergence of the NMR shielding 
with respect to the kinctic-cncrgy cutoff i?kin and the 
results are presented in Table I. The shielding a is con- 
verged to within 0.02 ppm for a kinetic-energy cutoff of 
80 Ryd. Similar tests on other structures show similar 
results. 



Next, we tested the convergence of the NMR shield- 
ing with respect to the magnitude of the microscopic 
dipole |ms| used and the energy convergence criterion. 
At first sight it might appear difhcult to accurately con- 
verge the electronic structure in the presence of a small 
microscopic magnetic point dipole. Thus, we tested using 
different magnitudes for the microscopic dipole spanning 
several orders of magnitude from 10~^ fiB (which is ac- 
tually much less than the value of a core spin) to 10^ 
(which is obviously much more than an electron spin). 
On the other hand, the ability to converge the electronic 
structure accurately goes hand in hand with the energy 
convergence criterion Ecom that is used in such calcula- 
tions. This criterion is defined such that the calculation 
is considered converged if the energy difference between 
two consecutive SCF steps is smaller than E'conv The 
results for the shielding as a function of jm,,! and Ecom 
are collected in Table II. It is interesting to see that it 
is just as simple to converge with a small dipole than 
it is to converge with a large dipole. In either case, us- 
ing at least E'conv = 10"** Ryd yields results converged 
to within 0.1 ppm. Such a convergence criterion is not 
even particularly "tight" and most standard codes use 
at least i?conv — 10^® Ryd as default. Note that first 
signs of non-linear effects appear if large dipoles such as 
|ms| = 100 /Ub or |ms| = 1000 /is are used. In conclu- 
sion of the above tests, we use |ms| = 1 /xb, iJconv = 10~® 
Ryd, and Skin = 100 Ryd for all calculations. 



A. Generation of GIPAW pseudopotentials 

We have generated special-purpose norm-conserving 
pseudopotentials for our GIPAW calculations. In addi- 
tion to standard norm-conserving pseudopotentials (PS) , 
the GIPAW pseudopotentials include (i) the full set of AE 
core atomic functions and (ii) the AE {(pn) and the PS 
(0„) valence atomic orbitals. The core orbitals contribu- 
tion to the isotropic NMR shielding is: 



1 E 2(2/„ + l) 



nGcore 



(27) 



The AE and PS valence orbitals arc used to compute 
the coefficients kji^nm and e-R,^nm at the beginning of the 
calculation. The PS valence orbitals are also used to 
compute the GIPAW projectors |pr) from: 



|PR,n) = '^{S ^)nm \4>'R,m) 
m 

Snm = {4>'R,n\4>Il,m) n 



(28) 
(29) 



S is the overlap between atomic PS wave function, inte- 
grated up to the cutoff radius of the corresponding pseu- 
dopotential channel. 

We construct at least two projectors per angular mo- 
mentum channel by combining each valence orbital with 
one excited state with the same angular momentum. For 
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TABLE II: Chemical shielding in ppm of a hydrogen atom in a water molecule. The shielding is given as a function of the 
magnitude of the microscopic dipole liiial (in units of Bohr magneton (xb) and the energy convergence criterion Econv (in units 
of Rydberg). -Econv is defined such that the calculation is considered converged if the energy difference between two consecutive 
SCF steps is smaller than Econv 





10"^ 


10"^ 


10"* 


10"^ 


10"® 


10-^ 


10"® 


io-« 


10-10 


0.00001 


31.5170 


31.2541 


31.2541 


31.1354 


31.1403 


31.1338 


31.1356 


31.1356 


31.1356 


0.0001 


31.5265 


31.2664 


31.2664 


31.1392 


31.1390 


31.1322 


31.1300 


31.1300 


31.1300 


0.001 


31.5253 


31.2667 


31.2667 


31.1395 


31.1397 


31.1328 


31.1304 


31.1304 


31.1304 


0.01 


31.5251 


31.2665 


31.2665 


31.1395 


31.1394 


31.1325 


31.1303 


31.1303 


31.1303 


0.1 


31.5250 


31.2664 


31.2664 


31.1393 


31.1393 


31.1323 


31.1301 


31.1301 


31.1301 


1.0 


31.5250 


31.2664 


31.2664 


31.1393 


31.1393 


31.1323 


31.1301 


31.1301 


31.1301 


10.0 


31.5250 


31.2663 


31.2663 


31.1392 


31.1392 


31.1322 


31.1301 


31.1301 


31.1301 


100.0 


31.5212 


31.2586 


31.2586 


31.1350 


31.1327 


31.1243 


31.1256 


31.1256 


31.1252 


1000.0 


31.0167 


30.5904 


30.7197 


30.6618 


30.6334 


30.6408 


30.6403 


30.6404 


30.6405 



TABLE III: Electronic configuration and cutoff radii for the 
norm-conserving pscudopotentials used in the present work. 



Atom 


configuration 


rc(s) 


rc{p) 


rc(d) 


H 


Is' 


0.50 






B 


[He] 2s' 2p^ 


1.40 


1.40 




C 


[He] 2s' 2p' 


1.50 


1.50 




N 


[He] 2s'2p^ 


1.45 


1.45 




O 


[He] 23^2/ 


1.40 


1.40 




F 


[He] 2s' 2p^ 


1.30 


1.30 




P 


[Ne] 3s'3p' ''3d'^ 


1.90 


2.10 


2.10 


Si 


[Ne] 3s'3p^-^3d° ' 


2.00 


2.00 


2.00 


CI 


[No] 3si-"''3/-53d°-"5 


1.40 


1.40 


1.40 


Cu 


[Ar] 4sUp°3d^° 


2.05 


2.20 


2.05 



IV. RESULTS 

In this section we present results for molecules and 
solids. We first calculated the absolute shielding tensor 
of some small molecules by two different approaches, the 
direct (linear response) and the converse method, in or- 
der to check that the two yield the same results. Then, 
we compared the chemical shifts of fluorine compounds, 
calculated by the converse method and by all-clcctron 
large basis set quantum-chemistry calculations. Finally, 
we report the calculated ^^Si chemical shifts of three Si02 
polymorphs and the Cu shift of a metallorganic com- 
pound. 

A. Small molecules 



example, for hydrogen we include the 2s orbital in the set 

of atomic wave functions. For all second row elements, 
we add the 3s and 3p orbitals, and so on. If any excited 
state turns out to be unbound (as in the case of oxygen 
and fluorine), we generate an atomic wave function as a 
scattering state at an energy 0.5 Ry higher than the cor- 
responding valence state. This procedure ensures that 
the GIPAW projectors are linearly independent and that 
the matrix S is not singular. 

We found that the accuracy of the calculated NMR 
chemical shifts depends critically on the cutoff radii of 
the pscudopotentials. Whereas the total energy and the 
molecular geometry converge more quickly with respect 
to reducing the pseudopotential radii, the NMR chem- 
ical shift converges more slowely. Therefore, GIPAW 
pscudopotentials have to be generated with smaller radii 
compared to the pscudopotentials usually employed for 
total energy calculations. Table III reports the atomic 
conflguration and the cutoff radii used to generate the 
pscudopotentials. 



We calculated the chemical shift of hydrogen, car- 
bon, fluorine, phosphorus, and silicon atoms of various 
small molecules. First, the structures were relaxed us- 
ing PWSCF in a box of 30 Bohr and a force convergence 
threshold of 10~^ Ry/Bohr. Using the relaxed positions, 
the chemical shifts were calculated using both the direct 
and converse method, and the results are shown in Ta- 
ble IV. This benchmark calculation shows that the direct 
and the converse methods agree to within less than 1%. 

B. Fluorine compounds 

In structural biology ^^F NMR spectroscopy plays an 
important role in determining the structure of protein 
membranes. The advantage over ^^N and ^^O label- 
ing is twofold: the natural abundance of "'^^F is nearly 
100%, and ^^F has spin 1/2, i.e. a vanishing nuclear 
quadrupole moment. Quadrupole interactions in high- 
spin nuclei (e.g. "'^^O) are responsible for the broadening 
of the NMR spectrum. On the contrary, ^^F NMR yields 
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TABLE IV: Results for the absoule NMR shielding a in ppm 
of small molecules for the direct and converse methods. The 
core contribution to the shielding is also shown. 



molecule 



H shielding 

CH4 

SiH4 
TMS 

C shielding 

CH4 
CeHe 
CH3F 
TMS 

F shielding 

CH3F 

PF3 

SiF4 

SiHgF 

P shielding 

P2 
PF3 

Si shielding 

SiF4 
SiH3F 

Si2H4 

TMS 



direct 



30.743 
22.439 
27.444 
30.117 

185.435 
36.887 
93.704 
175.774 

448.562 
277.148 
378.857 
423.456 

-323.566 
150.603 

431.438 
337.648 
230.830 
320.958 



30.670 
22.403 
27.413 
30.125 

186.027 
37.205 
94.250 

176.094 

447.014 
275.819 
376.227 
422.253 

-320.201 
150.856 

432.495 
337.677 
230.489 
320.636 



0.0 
0.0 
0.0 
0.0 

200.333 
200.333 
200.333 
200.333 

305.815 
305.815 
305.815 
305.815 

908.854 
908.854 

837.913 
837.913 
837.913 
837.913 
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FIG. 1: (Color online) Experimental vs. calculated ^®F chem- 
ical shifts in ppm, with respect to CF3CI. For sake of clarity 
the chemical shift of F2 is not shown. 



basis functions for 2nd row atoms. The calculation time 
required by our plane- wave converse method is compara- 
ble to that of cc-pV5Z calculations. 



very sharp and resolved lines. In addition, it has been 
found that the substitution of — CH3 groups with — CF3 
in some amino-acids does not perturb the structure and 
the activity of protein membranes, allowing for in vivo 
NMR measurements. 

In order to benchmark the accuracy of our method, 
we calculated the ^^F chemical shifts of ten fluorine com- 
pounds utilizing the converse method, and we compared 
our results to all-electron gaussian-basis set calculations, 
as well as to experimental data. The molecules were first 
relaxed with GaussianOS,^"* with the 6-311-|-g(2d,p) basis 
set at the B3LYP level. Then, we calculated the IGAIM 
chemical shift with GaussianOS, with the cc-pVTZ, cc- 
pVQZ, cc-pV5Z and cc-pV6Z basis sets.^^ 

To calculate the relative chemical shifts, we used CgFg 
as a secondary reference compound, and we used the ex- 
perimental CgFg chemical shift, to get the primary ref- 
erence absolute shift (CF3CI). The results are shown 
in Table V and in Fig. 1. While compiling Table V we 
suspected that the experimental values of P-C6H4F2 and 
CgHsF have been mistakenly exchanged in Ref. [26]. A 
quick inspection of the original paper, confirmed our 
suspicion. The overall agreement of the converse method 
with experimental data is very good, and of the same 
quality as the cc-pV6Z basis set, which comprises 140 



C. Solids 

In this section we present the ^^Si and "'^''O chem- 
ical shifts calculated by our converse method in four 
Si02 polymorphs: quartz, /?-cristobalite, coesite and 
stishovite. Coesite and stishovite are metastable phases 
that form at high temperature and pressures developed 
during a meteor impact. Besides their natural occur- 
rence in meteors, they can also be artificially synthesized 
by shock experiments. 

We adopted the experimental crystal structures and 
atom positions in all calculations. We used a cutoff of 
100 Ryd and a k-point mesh of 8x8x8 for quartz, /3- 
cristobalite and stishovite. In the case of coesite, having 
the largest primitve cell (48 atoms), we used a k-point 
mesh of 4x4x4. In Table VI we show a comparison be- 
tween the calculated and the experimental chemical shifts 
for the four crystals. We determined the ^^Si and ^^O 
reference shielding as the intercept of the least-square 
linear interpolation of the (cTcaic, <^cxpt) pairs. Note that 
the nuclear magnetic dipole rris breaks the symmetry of 
the hamiltonian. Thus, we retained only the symmetry 
operations that map site s in s' without changing the ori- 
entation of the magnetic dipole (i.e. s s' , nis = mg'). 

Another important point is that in periodic systems 
we are not just including one nuclear dipole, but rather 
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TABLE V: Experimental and calculated chemical shifts in ppm, with respect to CF3CI. In all calculations (GaussianOS 

and GIPAW plane waves), we used CfiFe as the reference compound. The exj)erimcntal values of P-C6H4F2 and CeHsF are 
exchanged with respect to Ref. [26]. For molecules with inequivalent F atoms. tlK> average chemical shift is reported. MAE is 
th(^ mean al)sohit(' error iu jipui. with resiiect to oxporimeiit. 



A/Tnlppii 1p 


Expt.^'^ 




(-raiicit;iari CC-T\\f 07i 

VJCtUOOlClll VjVj jj V Vi^ZJ 


(-rmiQciifin pp-ri\/Pi7i 








—251 






-254.25 


-9^4 70 


_oc;s ■?! 

^00.0 L 






—1 fi4 QO 


—1 fi4 Qfl 


-1 fi4 Qfl 


-1 fi4 Qfl 


—1 fi4 QO 


tits 


— ioi.o 


— i4o.9o 


— io9.oo 


—Lob.o 1 


1 /in 
-loo. 49 


1 Q K tic 


P-C6H4F2 


-113.15 


-115.77 


-113.98 


-114.07 


-113.78 


-111.84 


CeHsF 


-106 


-106.84 


-104.94 


-104.83 


-104.35 


-104.68 


(CF3)2CO 


-84.6 


-90.37 


-82.12 


-78.81 


-77.83 


-76.63 


CF3COOH 


-76.55 


-87.82 


-81.38 


-77.88 


-76.80 


-75.90 


CeHsCFa 


-63.72 


-78.24 


-69.42 


-66.21 


-65.28 


-64.16 


CF4 


-62.5 


-74.48 


-73.94 


-68.76 


-66.66 


-66.05 


F2 


422.92 


367.92 


375.36 


383.03 


385.72 


390.02 


MAE 




13.19 


9.40 


7.35 


6.69 


5.64 



TABLE VI: Calculated and experimental ^^Si and ^''O NMR 
chemical shifts of four Si02 crystals in ppm. Experimental 
data was taken from Refs. [28] and [29] . In the case of coesite 
all inequivalent chemical shifts are reported. 



Mineral 



quartz 

/3-cristobalite 

stishovite 
coesite 



quartz 

/3-cristobalite 

stishovite 

coesite 



Calc. 5 [ppm] 



Expt. 5 [ppm] 



^Si 



-107.10 
-108.78 

-184.13 
-107.30 
-113.35 

43.52 
40.35 
116.35 
26.35 
39.66 
52.70 
56.84 
59.03 



-107.73 
-108.50 

-191.33 
-107.73 
-114.33 

40.8 
37.2 

N/A 
29 
41 
53 
57 
58 



an infinite array. Thus, interactions between iiis and 
its infinite periodic replicas become important, and the 
chemical shift should be converged with respect to the 
superccU size. To test for this convergence, we repeated 
the calculations for quartz and /3-cristobalite in a larger 
supercell and we found a change in chemical shift of less 
than 0.1 ppm. This rapid convergence is due to the 
decay of the magnetic dipole interactions. 



D. Large Systems 



Reactive sites in biological systems such as 
organometallic molecules, as well as inorganic ma- 
terials, are of great importance. In particular, there is a 
surge of interest in studying copper(I) reactive sites using 
solid-state NMR. NMR experiments on these materials 
are challenging because of the large nuclear quadrupole 
moments of ^"^Cu and ^^Cu. Here, wc present the 
results for the copper-phosphine metallocene, tetram- 
ethylcyclopentadienyl copper(I) triphenylphosphine 
(CpCuPPhS), which as a solid contains 228 atoms in a 
primitive orthorhombic unit cell. The molecular struc- 
ture is shown in Fig. 2. The properties of the shielding 
tensor for the copper environment were observed exper- 
imentally for the solid material, and simulated using 
quantum-chemical methods on the molecular complex.^*^ 

While the converse approach can calculate the chemi- 
cal shift for this large system easily, it is more challeng- 
ing for the linear-response method, which in our experi- 
ence took much longer in general, did not finish at all, or 
was unable to handle such large systems. We calculated 
the copper chemical shift for CpCuPPhS using the con- 
verse method with an energy cutoff of 80 Ry in the self- 
consistent step and PBE pseudopotentials. While previ- 
ous quantum-chemical calculations were able to reason- 
ably reproduce the experimental span (1300 ppm) and 
the skew (0.95) of the chemical shielding tensor, they 
were not able to calculate the chemical shift itself (0 ppm 
relative to copper (I) chloride), with an inaccuracy of 
several hundred ppm.'^'^ In addition to yielding excellent 
agreement with experiment to within 2 ppm for the chem- 
ical shift, our calculations also gave good results for the 
span (1038 ppm) and the skew (0.82) of the chemical 
shielding tensor. 
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NMR method is built on Mauri's GIPAW method, but 
has the advantage of circumventing the need for a hnear 
response framework. 

The main advantage of our converse-NMR method is 
that it requires only the ground state wave functions 
and Hamiltonian to calculate the orbital magnetization. 
Since no external magnetic field is included in the cal- 
culation, our method solves the gauge-origin problem. 
Moreover, "difficult cases" can be treated easily by our 
converse method, provided that relativistic corrections 
and many-body effects are included in the Hamiltonian. 
Thus, one can concentrate effectively all efforts in de- 
veloping advanced post-DFT theories (i.e. DFT-I-U, 
DMFT, hybrid functionals, self-interaction-free methods) 
and benchmark them against NMR experiments. 



FIG. 2: (Color online) The molecular structure of the metal- 
locene, tetramethylcyclopentadienyl copper(I) triphenylphos- 
phine (CpCuPPhS). The copper is shown in gold (light gray) 
forming the metallocene bond, while the phosphorus is green 
(gray). The crystal structure of this materials consists of a 
228-atom orthorhombic unit cell. 



V. DISCUSSION 

The results presented in the previous sections and in 
a previous work^ show that DFT is able to predict accu- 
rately the chemical shift of molecules and solids. In gen- 
eral, we expect this to be true for any weakly-correlated 
system, well described by the generalized-gradient ap- 
proximation (GGA). In addition to that, relativistic cor- 
rections to the NMR chemical shifts are negligible for all 
light elements in the periodic table, and become impor- 
tant starting from 4th row elements. 

However, there will always exist "difficult" cases in 
which relativistic corrections cannot be neglected and/or 
one has to go beyond DFT with standard local function- 
als. This is an active field of research in quantum chem- 
istry^-^''^^ and today it is customary to compute NMR 
chemical shifts with semi-local hybrid DFT functionals 
(such as B3LYP). Most quantum-chemistry codes allow 
the inclusion of relativistic effects (spin orbit) by pertur- 
bation theory; furthermore, fuUy-relativistic (four com- 
ponent) solutions of the Dirac-Breit equation have re- 
cently been implemented.'^^ 

To the best of our knowledge, all existing ab-initio 
codes, calculate NMR shifts by perturbation theory. 
Among them, localized-basis sets are the most popu- 
lar choice to expand wave functions. This leads to 
very complicated mathematical expressions and to gauge- 
dependent results. Only two plane-wave, linear-response 
implementations^'^ have been reported. Our converse- 



VI. SUMMARY 

In this paper we have generalized the recently devel- 
oped converse NMR approach^ such that it can be used in 
conjunction with norm-conserving, non-local pseudopo- 
tentials. We have tested our approach both in finite and 
periodic systems, on small molecules, four silicate miner- 
als and a molecular crystal. In all cases, we have found 
very good agreement with established methods and ex- 
perimental results. 

The main advantage of the converse-NMR method is 
that it requires only the ground state wave functions and 
Hamiltonian, circumventing the need of any linear re- 
sponse treatment. This is of paramount importance for 
the rapid development and validation of new methods 
that go beyond DFT. 

Currently, we are applying the converse NMR method 
to study large biological systems such as nuclei acids'^^ 
and drug-DNA interactions'^^''^^ in conjunction with a 
recently developed van der Waals exchange-correlation 
functional. '^^''^^ We are also exploring the possibility to 
calculate non-perturbatively the Knight shift in metals. 
Finally, the converse method can be used to calculate the 
EPR g-tensor in molecules and solids.'^* 
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Appendix A: The GIPAW transformation 

The starting point is the projector augmented wave 
(PAW) transformation:^^ 



r=l + ^(|<^R,„)-|0R,„)) {m,n\ 



(Al) 



which connects an all-electron wave function |^) to the 

corresponding pseudopotential wave function l-tp) via: 
Itjj) ~ TIV")- Here, |(/>R,n) are all-electron partial waves, 
0R.ri) are pseudopotential partial waves, and (pr,„| are 
PAW projectors. The sum runs over the atom positions 
R. n is a combined index that runs over the set of projec- 
tors attached to atom R. In the original PAW formalism, 
there are two sets of projectors per angiilar momentum 
channel (0, . . . , ^max), each with {21 + 1) projectors for a 
total of 2(/inax + 1)^ PAW projectors. 

The expectation value of an all-electron operator Oae 
between all-electron wave functions can then be ex- 
pressed as the expectation value of a pseudo operator 
Ops between pseudo wave functions as (V'|Oae|V') = 
{i^\Ops\i^)j where the Ops is given by: 

Ops = T+OaeT = 

= Oae+ ^ \PR,n) {(pn,n\OA-E\<hl,m) - 
R,nm 

- (^R,n|OAE|^R,m) (PR,m| (A2) 

In the presence of external magnetic fields the PAW 
transformation is no longer invariant with respect to 
translations (except in the very simple case of only one 
augmentation region). This deficiency was resolved by 
Mauri et al. who developed the GIPAW ( "gauge includ- 
ing PAW") method^, which is similar to the PAW trans- 
formation from Eq. (A2) but with the inclusion of phase 
factor compensating the gauge term arising from the 
translation of a wave function in a magnetic field. The 
GIPAW transformation in the symmetric gauge reads: 



In the following we will refer to the often occurring part 
6ae = e-(*/2'=)'-^^^OAEe(^/2c)r-RxB -^^g^ operator 

and denote it with a hat. The above expression allows 
us to calculate accurate expectation values of operators 
within a pseudopotential approach. In this work we 
carry out the derivation working in the symmetric gauge 
A(r) = (1/2)B X r. This is not an issue, since all physical 
quantities we are working with, are gauge invariant. 

One useful property of the GIPAW transformation is: 



-(i/2c)r-RxB 



P + -A(r) 



„(j/2c)r-RxB _ 



p + iA(r-R)]" (A5) 



Appendix B: Periodic systems 



Note that the set of equations (23) are well defined also 

in periodic boundary conditions. In fact, Mbaro can be 
calculated by the Modern Theory of the Orbital Magne- 



tization 
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as: 



Mbare = -^^^'^fi^nk) {dkUnk\ 

nk 



X CHk + enk - 2eF) |9kWnk) 



GIPAW 



(Bl) 



Tg = 1 + ^ e(^/2^)-^x« (|.^R,„) - \m) ■ 

R,n 

-(i/2c)r-RxB 



(PR,i 



and the corresponding pseudopotential operator Ops is: 
Ops = Tg+OabTg = Oae + ^ e(i/2c)r.RxB |~^^^ . 

R,rtm 

■ (-/-R, Je-(V2^)-«x^ Oae e(^/2^)-^xB|^^_^^ _ 



- (.^R,n|e-(^/2c)-RxB g(i/2c)r.RxB|^^^^^ 



-(i/2c)r-RxB 



(A4) 



where Hk is the GIPAW hamiltonian, e„k and |w„k) are 
its eigenvalues and eigenvectors. 

The position operator appearing in the other terms 
in Eq. (23) is well defined because the projectors I/Jr,^) 
and 1pr,„) are non- vanishing only inside an augmentation 
sphere, centered around atom R. 

The expression of Mnl, Mpaia and Mpara can be fur- 
ther manipulated in order to work with Bloch wave func- 
tions. In the following part, we show it only for Mnl, 
since Mpara is similar. The term Mjia can be manipu- 
lated instead in a trivial way. 
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Let's consider the expectation value of 

(r - R) X [r - R, F^L] = _ r) ^ (yNL^ _ 

on a Bloch state IV'nk) = e*'*'' |M„k): 

- Yl Me-'^^'ir - R) X (r - R)e^''"-|^^nk) = (B3) 

R 

= - E E ("»k|e-*--(r - L - t) |/3L+.,i) x ^;,^.+- | (r - L - r)e*n«nk> , (B4) 

where L arc the real space lattice vectors (not to be confused with the angular momentum operator) and r are the 
position of the atoms in the unit cell. Inserting two canceling phase factors; 

• ■ • = - E E ("nk|e-*('-^-")(r - L - r) |/3L+x,i) x v^+^ (/^l+xjI (r - L - r)e*(--L-) |«„k) (B5) 

Lr ij 

one can recognize immediately the k derivative of the KB projectors. In addition, since the KB projectors vanish 
outside their augmentation regions, it is possible to insert a second sum over L' running on the right hand side of the 
cross product: 

■ • ■ = - E E («nk|(-l/i)ak (e-*(--^-) I^L+x,^)) X v^+^ (l/i)5k ((^L'+x,il e'''^'-^-")) M (B6) 

LL't ij 



r \ L / \ L' 



''(■■-^'-"^ |n„k) (B7) 



In periodic systems the structure factors can be absorbed by the projectors: 

< = J2J2\l^^,i)''r,ij{l3^,j\ (B8) 

T ij 

= l^e-''(-^--)|;3L+^,> (B9) 

L 

Finally: 



Mnl = 



^ E E 7 (^"k I 9k/3j^,i) X vr,ij (^k^J^,,- I u„k) (BlOa) 

nk T,ij 
^ occ ^ 

>^ E E 7 (""k I 9k^,j) X kr,ij (dkP^j I «„k) (BlOb) 

nk T,ij 
^ occ 

Mdia = -Y^{'^nk\p^^i)er,ij {p^j\Unl,) (BlOc) 



M 

para 



nk 



r 



This completes the main result and it allows us to calcu- efficiently calculate the NMR chemical shift for elements 
late the orbital magnetization in the presence of non-local heavier then hydrogen using Eq. (1). 
pseudopotentials. With this result we can now easily and 
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